POP Steering parameterization

Eddy diffusivity implementation by Bates, Tulloch, Marshall, and Ferrari

POP Steering Parameterization

CAM Implementation of Mesoscale eddy diffusivity based on mixing theory

Here are my notes on the development of a parameterization for eddy diffusivity. The scheme is based on the paper by Michael Bates, Ross Tulloch, John Marshall, AND Raffaele Ferrari. I'm working with Gokhan Danabasoglu and Matt Long here at NCAR and John Marshall at MIT.

Gokhan's Implementation Notes on Diffusivity with Steering Level Suppression

12 November 2014

Starting with $$K=u_{rms}∗L_{mix}$$

the general form for diffusivity, $K$, is given by equation (6) of [Bates et al. (2014)] http://journals.ametsoc.org/doi/abs/10.1175/JPO-D-13-0130.1. Namely,

$$L_{mix} = {\Gamma * L_{eddy} \over (1 + b1 * |u_{mean} - c|^2 /u_{rms}^2 (z=0)} $$

Here,

  • $u_{rms}$ is the root-mean-square (rms) eddy velocity;

  • $L_{mix}$ is the mixing length;

  • $L_{eddy}$ is the eddy diameter (depth independent);

  • $u_{mean}$ is the mean zonal velocity (resolved);

  • $c$ is the zonal eddy phase speed (depth independent);

  • $\Gamma = 0.35$;

  • $b1 \sim 4$.

Eddy Length Scales:

Rossby deformation radius = $L_r = {c_r \over |f|}$

Equatorial Rossby deformation radius $= L_{req} = {\sqrt{c_r \over 2\beta}}$

Rhines scale $= L_{Rh} = {\sqrt{u_{rms} \over \beta}} \sim {\sigma_{vi} \over \beta}$

$c_r$ is the first baroclinic wave speed computed following equation (2.2) of Chelton et al. (1998) with $m=1$;

$f$ is the Coriolis parameter;

$\beta$ is the latitudinal variation of the Coriolis parameter; and

$\sigma_{vi}$ is the Eady growth rate given by

$\sigma_{vi} = {f \over \sqrt{R_i}}$

with $R_i$ the vertically integrated (over the 100 – 2000 m depth range) Richardson number.

So, any one these length scales could be used as an eddy length scale. An alternative is

$L_{eddy} = min (L_r, L_{req}, L_{Rh}).$

Eddy Velocity:

$u_{rms} = alpha*\sigma_{vi}*L_r$

where $\sigma_{vi}$ is the Eady growth rate based on local Richardson number and $\alpha$ is a scaling constant.

Zonal phase speed:

$c = - \beta * L_r^2$

References:

Bates, M., R. Tulloch, J. Marshall, and R. Ferrari, 2014: Rationalizing the spatial distribution of mesoscale eddy diffusivity in terms of mixing length theory. J. Phys. Oceanogr., 44, 1523-1540, doi: 10.1175/JPO-D-13-0130.1.

Chelton, D. B., R. A. deSzoeke, M. G. Schlax, K. E. Naggar, and N. Siwertz, 1998: Geographical variability of the first baroclinic Rossby radius of deformation. J. Phys. Oceanogr., 28, 433-460.

Tullock, R., J. Marshall, and K. S. Smith, 2009: Interpretation of the propagation of surface altimetric observations in terms of planetary waves and geostropic turbulence. J. Geophys. Res., 114, C02005, doi: 10.1029/2008JC005055.

Questions:

  1. In the 2D implementation, $u_{rms}$ and $u_{mean}$ specifications: upper-ocean vertically or integrated or at $z = 0$?

    A: $U_{rms}$ is not depth dependent; both $u_{rms}$ and $u_{mean}$ are for surface only

  2. In the 2D implementation, vertical profile will be specified by $N2(z)$?

    A: Yes, $N^2(z) \over N_{ref}(z)$ to be more precise

  3. Local $R_i$ use imbedded in sigma in $u_{rms}$ calculation?

    A: No, $u_{rms}$ is depth independent

  4. alpha = ?

    A: trial and error

  5. Cancellation of $f$’s in $u_{rms}$ calculation?

  6. Zonal phase speed equation correct? Both $\beta$ and $L_r$ will be positive, producing $c < 0$ always. This appears to be in contrast with Tullock et al. (2009).

    A: What is plotted in Bates is (U-c)

Parameterizing the Eddy Length Scale

$$K=u_{rms}∗{\Gamma * \color{red}{L_{eddy}} \over (1 + b1 * |u_{mean} - c|^2 /u_{rms}^2 (z=0)} $$

The three length scales under consideration:

  1. Rossby Deformation Radius $\quad\quad L_r = {c_r \over |f|}$
  2. Rossby Equatorial Radius $ \quad\quad= L_{req} = {\sqrt{c_r \over 2\beta}}$
  3. Rhine's Scale $ \quad\quad L_{Rh} = {\sqrt{u_{rms} \over \beta}} \sim {\sigma_{vi} \over \beta}$

The Rossby Deformation Radius $L_r$ depends on the the Baroclinic wave speed $c_r$ and the Coriolis force $f$. The CESM value of the first Baroclinic wave speed is derived as per eq 2.2 in Chelton 1998.

First Baroclinic Wave Speed $c_r$ Chelton vs CESM

In [1]:
# <!-- collapse=True --> 
%pylab inline
%run steering_setup ;
basedir='/scratch/jet'
###nc = Dataset(basedir+'/steering/g.e11.GIAF.T62_gx1v6.steer.007.pop.h.nday1.yr.0249.nc')
nc = Dataset(basedir+'/steering/g.e11.GIAF.T62_gx1v6.steer.007.pop.h.nday1.0249-01.nc')
omega=7.2921e-5
deg2rad=pi/180
f=2*omega*sin(lat*deg2rad)
Populating the interactive namespace from numpy and matplotlib

/loan/jet/install/anaconda/lib/python2.7/site-packages/mpl_toolkits/__init__.py:2: UserWarning: Module argparse was already imported from /loan/jet/install/anaconda/lib/python2.7/argparse.pyc, but /loan/jet/install/anaconda/lib/python2.7/site-packages is being added to sys.path
  __import__('pkg_resources').declare_namespace(__name__)

In [2]:
# <!-- collapse=True --> 
### get Variable from file c_rossby is in cm/s convert to m/s
c_rossby = nc.variables['C_ROSSBY'][0,:,:]


#Resample (aka re-project, re-grid) the NCEP data to target grid. First with nearest neighbour resampling...
c_rossby_nearest = pyresample.kd_tree.resample_nearest(orig_def, c_rossby, \
        targ_def, radius_of_influence=500000, fill_value=None)
c_rossby_mps = c_rossby_nearest/100.

fig1=MapZoneContour(lon_targcyc,lat_targcyc,c_rossby_mps,figsize=(14,9),
                    title="CESM First Baroclinic Wave Speed (m/s)",
                    fmt='%.1f' )
c_r_chelton=Image(filename=basedir+'/steering/chelton_sfig1.jpg')
display(c_r_chelton)
/loan/jet/install/anaconda/lib/python2.7/site-packages/matplotlib/text.py:52: UnicodeWarning: Unicode equal comparison failed to convert both arguments to Unicode - interpreting them as being unequal
  if rotation in ('horizontal', None):
/loan/jet/install/anaconda/lib/python2.7/site-packages/matplotlib/text.py:54: UnicodeWarning: Unicode equal comparison failed to convert both arguments to Unicode - interpreting them as being unequal
  elif rotation == 'vertical':

First Baroclinic Rossby Radius Chelton vs CESM

CESM Rossby deformation radius = $L_{eddy}$

$L_r = {c_r \over |f|}$

$L_{req} = {\sqrt{c_r \over 2\beta}}$

$\require{cancel} \cancel {L_{Rh} = {\sqrt{u_{rms} \over \beta}} \sim {\sigma \over \beta}}$

$L_{eddy} = min (L_r, L_{req}, \cancel{L_{Rh}}).$

In [3]:
# <!-- collapse=True --> 
fcor = nc.variables['FCOR'][0,:,:]
fcor_nearest = pyresample.kd_tree.resample_nearest(orig_def, fcor, \
        targ_def, radius_of_influence=500000, fill_value=None)
absf=np.abs(fcor_nearest)
l_rossby_calc=c_rossby_nearest/absf

btp = nc.variables['BTP'][0,:,:]
btp_nearest = pyresample.kd_tree.resample_nearest(orig_def, btp, \
        targ_def, radius_of_influence=500000, fill_value=None)

l_rossbyeq_calc=np.sqrt(c_rossby_nearest/(2*btp_nearest))
l_eddy_calc=np.minimum(l_rossby_calc,l_rossbyeq_calc)

### convert regridded radius from cm to km
cm2km=1.e-5
l_eddy_calc_km=l_eddy_calc*cm2km
clevs=[10,20,30,40,50,60,80,100,150,230]
fig2=MapZoneContour(lon_targcyc,lat_targcyc,l_eddy_calc_km,figsize=(14,9),
                    levels=clevs,
                    title="CESM First Baroclinic Rossby Radius (km) $L_{eddy} \sim min(l_r,l_{req})$",
                    fmt='%.1i' )
l_r_chelton=Image(filename=basedir+'/steering/chelton_fig2.jpg')
display(l_r_chelton)
/loan/jet/install/anaconda/lib/python2.7/site-packages/numpy/ma/core.py:937: RuntimeWarning: overflow encountered in multiply
  result = self.f(da, db, *args, **kwargs)

Parameterizing the Zonal Eddy Phase Speed $\color{red}{c}$

$K=u_{rms}∗{\Gamma * L_{eddy} \over (1 + b1 * |u_{mean} - \color{red}{c}|^2 /u_{rms}^2 (z=0)} $

$\cancel{c = - \beta * {L_r^2}}$ $L_r$ too high at equator

$c = - \beta * L_{eddy}^2$

In [4]:
# <!-- collapse=True --> 
#beta
clevs=arange(1.e-14,25.e-14,1e-14)
fig,ax,cbar =MapContourg(lon_targcyc,lat_targcyc,btp_nearest,
                    levels1=clevs,
#                    norm=norm,
                    addzonal=True,
                    figsize=(14,5),
                    title="Beta (1/cm*s)")
l_eddy_calc_sq=l_eddy_calc*l_eddy_calc
clevs=np.logspace(1,15,40)
norm=LogNorm()
fig,ax,cbar =MapContourg(lon_targcyc,lat_targcyc,l_eddy_calc_sq,
                    addzonal=True,
                    levels1=clevs,
                    norm=norm,
                    figsize=(14,5),
                    title="$L_{eddy}^2 (cm^2)$ log scale")
#cbar.ax.get_xaxis().get_major_formatter().labelOnlyBase = True
#cbar.ax.get_yaxis().set_major_formatter(plt.LogFormatter(10,  labelOnlyBase=True))
#ax.yaxis.set_major_locator(ticker.MultipleLocator)
#ax.yaxis.set_major_formatter(LogFormatter())
#cbar.formatter.set_scientific(False)
#cbar.formatter.set_powerlimits((0, 0))
#cbar.update_ticks()
cbar.set_ticks(np.logspace(1,15,15))

c=btp_nearest*l_eddy_calc*l_eddy_calc # units=1/(cm s) * cm^2 = cm/s


clevs=arange(-10,20,1)
fig,ax,cbar = MapContourg(lon_targcyc,lat_targcyc,c,
                    addzonal=True,
                    levels1=clevs,
#                    norm=norm,
                    setover='darkred',
                    figsize=(14,5),
                    title="eddy phase speed c = beta*leddy**2 (cm/s) linear scale(NOTE I'm plotting positive c here - changed to negative after this plot")
cbar.set_ticks(np.linspace(-10,20,7))
c=-1.*c
/loan/jet/install/anaconda/lib/python2.7/site-packages/matplotlib/contour.py:1515: UserWarning: Log scale: values of z <= 0 have been masked
  warnings.warn('Log scale: values of z <= 0 have been masked')

Hughes Phase Speed (cm/s) from Tulloch Marshall Smith '09

In [5]:
# <!-- collapse=True --> 
hughes_c=Image(filename=basedir+'/steering/Hughes_phase_speed.png')
display(hughes_c)

Parameterizing the Zonal Velocity $\color{red}{u_{mean}}$

$K=u_{rms}∗{\Gamma * L_{eddy} \over (1 + b1 * |\color{red}{u_{mean}} - c|^2 /u_{rms}^2 (z=0)} $

In [6]:
# <!-- collapse=True --> 
#################  PLOT (u-c)^2  convert m^2/s^2 #######################
umean=nc.variables['U_MEAN'][0,:,:]
umean10_nearest = pyresample.kd_tree.resample_nearest(orig_def, umean,
                targ_def, radius_of_influence=500000, fill_value=None)

umean10_mps=umean10_nearest*1.e-2

clevs=np.linspace(-60,60,40)
fig,ax,cbar=MapContourg(lon_targcyc,lat_targcyc,umean10_nearest,
                    addzonal=True,
                    figsize=(14,5),
                    levels1=clevs,
#                    setover='darkred',
#                    setunder='darkblue',
                    title="U zonal velocity (cm/s)")
#fig.subplots_adjust(left=0.15, right=.9, bottom=0.1, top=0.9);
cbar.set_ticks(arange(-60,60,10))
pylab.xlim([-20.,20.])


# get 2D versions of the lat and lon variables add longitude start point here!
lon_ecco_orig, lat_ecco_orig = np.meshgrid(lon_ecco[:], lat_ecco[:])

orig_ecco_def = pyresample.geometry.GridDefinition(lons=lon_ecco_orig, lats=lat_ecco_orig)

ecco_u_nearest = pyresample.kd_tree.resample_nearest(orig_ecco_def, ecco_u, \
        targ_def, radius_of_influence=500000, fill_value=None)
ecco_u_nearest_masked=np.ma.masked_where(ecco_u_nearest<-10, ecco_u_nearest, copy=True)

#print ecco_u_cyc_masked[:,0]

clevs=np.linspace(-60,60,40)
#ecco_u_nearest = pyresample.kd_tree.resample_nearest(orig_def, ecco_u, \
#        targ_def, radius_of_influence=500000, fill_value=None)
fig,ax,cbar=MapContourg(lon_targcyc,lat_targcyc,ecco_u_nearest_masked*100,
                    addzonal=True,
                    figsize=(14,5),
                    levels1=clevs,
#                    extend='both',
                    title="ECCO U zonal velocity (cm/s)")
cbar.set_ticks(arange(-60,60,10))
pylab.xlim([-20.,20.])
Out[6]:
(-20.0, 20.0)
In [7]:
# <!-- collapse=True --> 
clevs=-np.logspace(3,-3,20)
norm=matplotlib.colors.SymLogNorm(linthresh=1e-4, vmin=-1000, vmax=-.001)
fig,ax,cbar=MapContourg(lon_targcyc,lat_targcyc,c,
                    addzonal=True,
                    levels1=clevs,
                    norm=norm,
                    figsize=(14,5),
                    title="eddy phase speed c = beta*leddy**2 (cm/s) log scale")

# Customize y tick lables
yticks=[-1000,-500,-200,-100,-50,-20,-10,-5,-2,-1,-.500,-.200,-.1,-.05,-.02,-.01,-.005,-.002,-.001]
cbar.set_ticks(yticks)
cbar.set_ticklabels(yticks)

#ax.ticklabel_format(axis='y', style='sci', scilimits=(-2,2))
#majorFormatter = ticker.FormatStrFormatter('%.3f')
#majorLocator   = ticker.LogLocator(base=10.0, subs=[1.0], numdecs=6, numticks=15)
#ticker.MultipleLocator(20)
#cbar.ax.yaxis.set_major_locator(majorLocator)
#cbar.ax.yaxis.set_major_formatter(majorFormatter)
#ax.yaxis.get_major_formatter().set_powerlimits((0, 1))
pylab.xlim([-200.,0.])



uminusc=umean10_nearest-c  #units cm/s
clevs=np.linspace(0,25,25)
fig,ax,cbar=MapContourg(lon_targcyc,lat_targcyc,uminusc,
                    addzonal=True,
#                    norm=norm,
                    levels1=clevs,
                    figsize=(14,5),
                    title="U-c (cm/s) (max = %i/min = %i) (c calculated using Leddy)"%(uminusc.max(),uminusc.min()))
yticks=linspace(0,25,6)
zonalticks=np.linspace(0,20,3)
cbar.set_ticks(yticks)
cbar.set_ticklabels(yticks)
pylab.xlim([0,25])
pylab.xticks(zonalticks)
Out[7]:
([<matplotlib.axis.XTick at 0x2b3e7ca01290>,
  <matplotlib.axis.XTick at 0x2b3e7ca010d0>,
  <matplotlib.axis.XTick at 0x2b3e7d0589d0>],
 <a list of 3 Text xticklabel objects>)
In [8]:
# <!-- collapse=True --> 
bates_uminusc=Image(filename=basedir+'/steering/bates2013_uminusc.png')
display(bates_uminusc)
In [9]:
# <!-- collapse=True --> 
uminusc_mps_sq=uminusc*uminusc*1.e-4   #units cm/s*cm/s*m/100cm*m/100cm
clevs=np.logspace(-4,1,40)
norm=LogNorm()
fig,ax,cbar=MapContourg(lon_targcyc,lat_targcyc,uminusc_mps_sq,
                    addzonal=True,
                    levels1=clevs,
                    figsize=(14,5),
                    norm=norm,
                    title=" ${(U-c)}^2 {(m^2/s^2)}$ max = %g/min = %g (log scale) "%(uminusc_mps_sq.max(),uminusc_mps_sq.min()))
#customize colorbar ticks
yticks=[1e-4,3e-4,1e-3,3e-3,1e-2,3e-2,1e-1,3e-1,1,3,10]
cbar.set_ticks(yticks)
#customize zonal ticks
pylab.xlim([0,.06])
pylab.xticks([0,.03,.06])

clevs=np.linspace(0,.06,25)
fig,ax,cbar=MapContourg(lon_targcyc,lat_targcyc,uminusc_mps_sq,
                    addzonal=True,
                    levels1=clevs,
                    figsize=(14,5),
                    title=" ${(U-c)}^2 {(m^2/s^2)}$ max = %g/min = %g "%(uminusc_mps_sq.max(),uminusc_mps_sq.min()))
#customize colorbar ticks
yticks=[0,.02,.04,.06]
cbar.set_ticks(yticks)
#customize zonal ticks
pylab.xlim([0,.06])
pylab.xticks([0,.03,.06])
Out[9]:
([<matplotlib.axis.XTick at 0x2b3e7c5ddf90>,
  <matplotlib.axis.XTick at 0x2b3e7c5d35d0>,
  <matplotlib.axis.XTick at 0x2b3e7c451490>],
 <a list of 3 Text xticklabel objects>)
In [10]:
# <!-- collapse=True --> 
bates_uminuscsq=Image(filename=basedir+'/steering/Bates_u_minus_c_sq.jpg')
display(bates_uminuscsq)

Parameterizing the root-mean-square eddy velocity $\color{red}{u_{rms}}$ and the Eady Growth Rate $\color{red}{\sigma_{vi}}$

$K=\color{red}{u_{rms}}∗{\Gamma * L_{eddy} \over (1 + b1 * |u_{mean} - c|^2 /u_{rms}^2 (z=0)} $

$u_{rms} = alpha*\color{red}{\sigma_{vi}}*L_{eddy}$

The original derivation of $\sigma_{vi}$:

$\sigma_{vi} = {f \over \sqrt{R_i}}$

with $R_i$ the vertically integrated (over the 100 – 2000 m depth range) Richardson number.

An alternate derivation of $\sigma_{vi}$ is:

$ R_i = {N^2 \over { ( \frac {\partial u}{\partial z} )^2 + ( \frac {\partial v}{\partial z} )^2} }$

$ N^2 = {{-g \over \rho_0 }\frac {\partial \rho}{\partial z}} $

After hydrostatic and geostrophic approximations

$f \frac {\partial v}{\partial z} = {{-g \over \rho_0 }\frac {\partial \rho}{\partial x}}; \quad\quad f \frac {\partial u}{\partial z} = {{g \over \rho_0 }\frac {\partial \rho}{\partial y}} $

so

$\frac {\partial v}{\partial z} = {{-1\over f}{g \over \rho_0 }\frac {\partial \rho}{\partial x}}; \quad\quad \frac {\partial u}{\partial z} = {{1\over f}{g \over \rho_0 }\frac {\partial \rho}{\partial y}} $

$\therefore$

$ R_i = {f^2 N^2 \over \underbrace{{ {g^2 \over \rho_0^2 } ( \frac{\partial \rho}{\partial y})^2 + {g^2 \over \rho_0^2 } ( \frac{\partial \rho}{\partial z})^2} }_\text{m^4}}\quad\quad = \quad\quad {f^2N^2 \over m^4} $

$\sigma_{vi} = {f \over \sqrt{R_i}}\quad = \quad {f \over \sqrt{{f^2N^2 \over m^4}}}\quad = \quad {{\cancel f m^2} \over \cancel f N}$

$RX_1 = RX_{east} = \Delta\rho_x = \rho_{i+1,j} - \rho_{i,j}$

$RY_1 = RY_{north} = \Delta\rho_y = \rho_{i,j+1} - \rho_{i,j}$

$RZ_1 = RZ_{k+1} = \Delta\rho_z = \rho_{k} - \rho_{k+1}$

$\displaystyle{1 \over L_{R_i}} \displaystyle\int_{2000m}^{100m} \left\lbrace { {-g\over\rho_0}{\frac {\partial \rho} {\partial z}} \over { {g^2 \over \rho_0^2 } \left[( \frac{\partial \rho}{\partial y})^2 + ( \frac{\partial \rho}{\partial z})^2\right] } } \right\rbrace dz$

Note: missing $f^2$ which will be cancelled when forming $\sigma_{vi}$

$\quad\quad$ so $\cdots$ this is not $R_i$

Implementation notes

For Level K

Numerator : Top $= -grav * RZ_{SAVE}(\cdots k+1) * dzwr(k)$

Denominator :

$ \begin{align} work1 = p25 & * ( RX(..,i_{east},k)^2 \\ & + RX(..,i_{west},k)^2 \\ & + RX(..,i_{east},k+1)^2 \\ & + RX(..,i_{west},k+1)^2 ) / DXT(i,j)^2 \\ \end{align} $

$ \begin{align} work2 = p25 & * ( RY(..,j_{north},k)^2 \\ & + RY(..,j_{south},k)^2 \\ & + RY(..,j_{north},k+1)^2 \\ & + RY(..,j_{south},k+1)^2 \\ \end{align} $

$work3 = {\left( TOP \over (grav^2*(work1+work2))\right)}*dzw(k)$

Notes:

1)Need to be careful at top and bottom of ocean
2)Accurate dzw(k) for each (i,j) to form $L_{R_i}$
3) When constructing $sigma$ itself, use $RZ_{SAVE}$ with a minimum N value
4) use eps2
In [11]:
# <!-- collapse=True --> 
sigma_old = nc.variables['SIGMA_AVG'][0,:,:]
sigma_old_nearest = pyresample.kd_tree.resample_nearest(orig_def, sigma_old, \
        targ_def, radius_of_influence=500000, fill_value=None)

sigma_avg1 = nc.variables['SIGMA_AVG1'][0,:,:]
sigma_avg1_nearest = pyresample.kd_tree.resample_nearest(orig_def, sigma_avg1, \
        targ_def, radius_of_influence=500000, fill_value=None)
clevs=arange(.005,.05,.0025)
fig,ax,cbar=MapContourg(lon_targcyc,lat_targcyc,sigma_avg1_nearest*86400,
                levels1=clevs,
                figsize=(14,5),
                setunder='darkblue',
                setover='darkred',
                title="New Eady inverse time scale calc using alt sigma $(days^{-1}) \sim  {{m^2} \over N}$")
yticks=[.01,.02,.03,.04,.05]
cbar.set_ticks(yticks)

tulloch_eady=Image(filename=basedir+'/steering/tulloch_eady.png')
display(tulloch_eady)
In [12]:
# <!-- collapse=True --> 
#################  PLOT u_rms with Leddy #######################
cm2m=.01
alpha=5.
u_rms_leddy = alpha*sigma_avg1_nearest*l_eddy_calc    # units = 1*1/s*cm = cm/s
u_rms_leddy_mps=u_rms_leddy*cm2m

clevs=np.linspace(0,.4,40)
fig,ax,cbar=MapContourg(lon_targcyc,lat_targcyc,u_rms_leddy_mps,
                    addzonal=True,
                    levels1=clevs,
                    figsize=(14,5),
                    title="$u_{rms} ({m \over s}) \quad alpha * \sigma * l_{eddy}$ (alpha=%i,max = %i/min = %i) "%(alpha,u_rms_leddy_mps.max(),u_rms_leddy_mps.min()))
#customize colorbar ticks
yticks=[0,.1,.2,.3,.4]
cbar.set_ticks(yticks)
#customize zonal ticks
pylab.xlim([0,.2])
pylab.xticks([0,.1,.2])
Out[12]:
([<matplotlib.axis.XTick at 0x2b3e7db7e150>,
  <matplotlib.axis.XTick at 0x2b3e7dcdd290>,
  <matplotlib.axis.XTick at 0x2b3e7d8f7a50>],
 <a list of 3 Text xticklabel objects>)
In [13]:
# <!-- collapse=True --> 
bates_urms=Image(filename=basedir+'/steering/bates2013-urms.png')
display(bates_urms)
In [14]:
# <!-- collapse=True --> 
#################  PLOT u_rms^2  convert m^2/s^2 #######################
u_rms_leddy_mps_sq=u_rms_leddy_mps*u_rms_leddy_mps
clevs=np.linspace(0,.06,40)
norm=LogNorm()
fig,ax,cbar=MapContourg(lon_targcyc,lat_targcyc,u_rms_leddy_mps_sq,
                    addzonal=True,
                    levels1=clevs,
                    figsize=(14,5),
                    title="$u_{rms}^2 ({m^2 \over s^2})$ (alpha=%i,max = %g/min = %g) "%(alpha,u_rms_leddy_mps_sq.max(),u_rms_leddy_mps_sq.min()))
#customize colorbar ticks
yticks=[0,.02,.04,.06]
cbar.set_ticks(yticks)
#customize zonal ticks
pylab.xlim([0,.06])
pylab.xticks([0,.03,.06])
Out[14]:
([<matplotlib.axis.XTick at 0x2b3e7e51a090>,
  <matplotlib.axis.XTick at 0x2b3e7e51a510>,
  <matplotlib.axis.XTick at 0x2b3e7e12bb50>],
 <a list of 3 Text xticklabel objects>)
In [15]:
# <!-- collapse=True --> 
urms_bates=Image(filename=basedir+'/steering/Bates_urms_sq.jpg')
display(urms_bates)
In [16]:
# <!-- collapse=True --> 
uminusc_sq_ovr_urms_sq = uminusc_mps_sq/u_rms_leddy_mps_sq             #units cm/s /  cm/s  
uminusc_sq_ovr_urms_sq_masked=np.ma.masked_where(uminusc_sq_ovr_urms_sq >100., uminusc_sq_ovr_urms_sq, copy=True)

clevs=np.logspace(-3,4,21)
norm=LogNorm()
fig,ax,cbar=MapContourg(lon_targcyc,lat_targcyc,uminusc_sq_ovr_urms_sq_masked,                
                    addzonal=True,
                    levels1=clevs,
                    norm=norm,
                    figsize=(14,5),
                    setover='darkred',
                    setunder='darkblue',
                    title="${(U-c)}^2 \over u_{rms}^2$ (masking values over 100)")                    
#customize colorbar ticks
yticks=np.logspace(-3,4,15)
cbar.set_ticks(yticks)
#customize zonal ticks
pylab.xlim([0,100])
pylab.xticks([0,50,100])

clevs=arange(0,5,.1)
fig,ax,cbar=MapContourg(lon_targcyc,lat_targcyc,uminusc_sq_ovr_urms_sq_masked,                
                    addzonal=True,
                    levels1=clevs,
                    figsize=(14,4),
                    setover='darkred',
                    setunder='darkblue',
                    title="${(U-c)}^2 \over u_{rms}^2$ (masking values over 100)")
#customize colorbar ticks
yticks=[0,2,4]
cbar.set_ticks(yticks)
#customize zonal ticks
pylab.xlim([0,8])
pylab.xticks([0,4,8])
Out[16]:
([<matplotlib.axis.XTick at 0x2b3e7f25d450>,
  <matplotlib.axis.XTick at 0x2b3e7f25d490>,
  <matplotlib.axis.XTick at 0x2b3e7efce8d0>],
 <a list of 3 Text xticklabel objects>)
In [17]:
# <!-- collapse=True --> 
bates_uminusc_sq_ovr_urms_sq =Image(filename=basedir+'/steering/Bates_suppress.jpg')
display(bates_uminusc_sq_ovr_urms_sq )

Parameterizing the Mixing term $\color{red}{L_{mix}}$

$\color{red}{L_{mix}} = \Gamma * L_{eddy} * Suppression$

$\Gamma = 0.35$

$b1 = 4$

$Suppression= {1 \over (1 + b1 * |u_{mean} - c|^2 /u_{rms (z=0)}^2 )}$

$\color{red}{L_{mix}} = {\Gamma * L_{eddy} \over (1 + b1 * |u_{mean} - c|^2 /u_{rms}^2 (z=0)}$

In [18]:
# <!-- collapse=True --> 
gamma=.35
b1=4.0
supp=1./(1+b1*uminusc_sq_ovr_urms_sq)
supp1b1=1./(1+1*uminusc_sq_ovr_urms_sq)

cm2m=.01
l_eddy_calc_m=l_eddy_calc*cm2m
lmix=gamma*l_eddy_calc_m*supp
lmix1b1=gamma*l_eddy_calc_m*supp1b1


clevs=np.logspace(0,5,40)
norm=LogNorm()
fig,ax,cbar=MapContourg(lon_targcyc,lat_targcyc,lmix,
                    addzonal=True,
                    norm=norm,
                    levels1=clevs,
                    figsize=(14,5),
                    title="$L_{mix} = {\Gamma * L_{eddy} * Suppression} $")

pylab.xlim([0,10e3])
cbar.set_ticks([1,3,1e1,3e1,1e2,3e2,1e3,3e3,1e4,3e4,1e5])



clevs=np.linspace(0,1,30)
norm=matplotlib.colors.Normalize()
fig,ax,cbar=MapContourg(lon_targcyc,lat_targcyc,supp,
                    addzonal=True,
                    levels1=clevs,
                    figsize=(14,5),
                    title="Suppression factor ${1 \over (1 + b1 * |u_{mean} - c|^2 /u_{rms (z=0)}^2)}$")

#customize colorbar ticks
yticks=np.linspace(0, 1, 6)
cbar.set_ticks(yticks)
#customize zonal ticks
pylab.xlim([0,.8])
pylab.xticks([0,.25,.5,.75])
Out[18]:
([<matplotlib.axis.XTick at 0x2b3e986ff510>,
  <matplotlib.axis.XTick at 0x2b3e7e4f5510>,
  <matplotlib.axis.XTick at 0x2b3e9831fe10>,
  <matplotlib.axis.XTick at 0x2b3e7f3d4cd0>],
 <a list of 4 Text xticklabel objects>)
In [19]:
# <!-- collapse=True --> 
bates_suppression =Image(filename=basedir+'/steering/bates2013_suppression.png')
display(bates_suppression )

Parameterizing Eddy diffusivity $\color{red}{K}$

$\color{red}{K}=u_{rms}∗L_{mix}$

(b1, the suppression scaling factor, defaults to 4 but I also included a $\color{red}{K}$ plot with b1=1)

In [20]:
# <!-- collapse=True --> 
kdiff=u_rms_leddy_mps*lmix
kdiff1b1=u_rms_leddy_mps*lmix1b1

clevs=np.logspace(0,4,40)
norm=LogNorm()
fig,ax,cbar=MapContourg(lon_targcyc,lat_targcyc,kdiff,
                    addzonal=True,
                    norm=norm,
                    levels1=clevs,
                    figsize=(14,5),
                    title="Eddy diffusivity $K$ (log scale)")
#customize colorbar ticks
yticks=np.logspace(0, 4,5)
cbar.set_ticks(yticks)
#customize zonal ticks
pylab.xlim([0,11e3])
pylab.xticks([0,5e3,10e3])

clevs=np.linspace(0,1.e4,40)
fig,ax,cbar=MapContourg(lon_targcyc,lat_targcyc,kdiff,
                    addzonal=True,
                    levels1=clevs,
                    figsize=(14,5),
                    title="Eddy diffusivity $K$ ")
#customize colorbar ticks
#cbar.set_ticks([1,3,1e1,3e1,1e2,3e2,1e3,3e3,1e4,])
yticks=np.linspace(0, 10000,6)
print yticks
cbar.set_ticks(yticks)
#customize zonal ticks
pylab.xlim([0,11e3])
pylab.xticks([0,5e3,10e3])


clevs=np.linspace(0,1.e4,40)
fig,ax,cbar=MapContourg(lon_targcyc,lat_targcyc,kdiff1b1,
                    addzonal=True,
                    levels1=clevs,
                    figsize=(14,5),
                    title="Eddy diffusivity $K$ With a suppression scaling factor $b1$ of 1 instead of 4")
#customize colorbar ticks
#cbar.set_ticks([1,3,1e1,3e1,1e2,3e2,1e3,3e3,1e4,])
yticks=np.linspace(0, 10000,6)
print yticks
cbar.set_ticks(yticks)
#customize zonal ticks
pylab.xlim([0,11e3])
pylab.xticks([0,5e3,10e3])
[     0.   2000.   4000.   6000.   8000.  10000.]
[     0.   2000.   4000.   6000.   8000.  10000.]

Out[20]:
([<matplotlib.axis.XTick at 0x2b3e9955ff10>,
  <matplotlib.axis.XTick at 0x2b3e9955fa10>,
  <matplotlib.axis.XTick at 0x2b3e9927ffd0>],
 <a list of 3 Text xticklabel objects>)
In [21]:
# <!-- collapse=True --> 
bates_k =Image(filename=basedir+'/steering/bates2013_K.png')
display(bates_k )

Steering Level status

1. L_eddy values ok
2. Rossby wave speed ok
3. Eady inverse timescale ok
    a. Use alternate calculation otherwise has problems at equator and
    b.Richardson averages are too high causing sigma to drop too low
4. Zonal Eddy Phase Speed (c) 
    a. Too high 10S:10N
    b. Missing westward velocities in ACC which appear 
        in Hughes Data used in this parameterization
5. (U-c) Because of problem with c this term is also too high 10S:10N.
6. U_rms  (alpha*sigma*l_eddy)  OK        
    a. Using alpha of 5
7. Zonal mean velocity looks ok compared to ECCO annual average
8. Suppression OK
9. K is low by a factor of 5

Todo

1. Need to work on Zonal phase speed c
2. Plot vertical distribution of K
In [22]:
# <!-- collapse=True --> 
urms_bates2=Image(filename=basedir+'/steering/Bates-Tulloch-etal-2014-002.jpg')
urms_bates3=Image(filename=basedir+'/steering/Bates-Tulloch-etal-2014-003.jpg')
urms_bates4=Image(filename=basedir+'/steering/Bates-Tulloch-etal-2014-004.jpg')
urms_bates5=Image(filename=basedir+'/steering/Bates-Tulloch-etal-2014-005.jpg')
urms_bates6=Image(filename=basedir+'/steering/Bates-Tulloch-etal-2014-006.jpg')
urms_bates7=Image(filename=basedir+'/steering/Bates-Tulloch-etal-2014-007.jpg')
urms_bates8=Image(filename=basedir+'/steering/Bates-Tulloch-etal-2014-008.jpg')
urms_bates9=Image(filename=basedir+'/steering/Bates-Tulloch-etal-2014-009.jpg')
urms_bates10=Image(filename=basedir+'/steering/Bates-Tulloch-etal-2014-010.jpg')
urms_bates11=Image(filename=basedir+'/steering/Bates-Tulloch-etal-2014-011.jpg')
display(urms_bates3,urms_bates4,urms_bates5,urms_bates6,urms_bates7,urms_bates8,urms_bates9,urms_bates10,urms_bates11)

CESM $U_{mean} = U_{resolved}$ integraged over {10,50,100,200,500} meters

In [23]:
# <!-- collapse=True --> 
umean_cesm=Image(filename=basedir+'/steering/umean.jpg')
display(umean_cesm)
In [24]:
# <!-- collapse=True --> 
bates_len_vel=Image(filename=basedir+'/steering/bates_len_vel.jpg')
display(bates_len_vel)